DALS028-批次效应04-因子分析

发表于 | 分类于
| 字数统计: 3,761 | 阅读时长 ≈ 17

前言

这一部分是《Data Analysis for the life sciences》的第10章批次效应的第3小节，这一部分的主要内容涉及批次效应(Batch Effects)的校正，有关这一部分Rmarkdown文档参见作者的Github。

因子分析

在我们介绍另一种用于批量效应校正的统计方法之前，我们先介绍这种方法的核心统计思想：因子分析(Factor Analysis)。因子分子早是在一个多世纪前出现的。卡尔·皮尔森(Karl Pearson)指出，当计算学生之间不同科目的数据时，这些不同科目之间存在着相关性。为了解释这种现象，他提出了一个模型，该模型有一个能解释这种相关性的因子，对于每个学生来说，这个因子在各个学科中都是通用的(common)，如下所示：

其中，$Y_{ij}$ 表示每个人 $i$ 的科目 $j$ 的成绩，而$\alpha_{i} $表示学生 $i$ 获得好成绩的能力。

在这个案例中，$W_{1}$ 是一个常数。这里我们使用一个稍微复杂的情况来说明因子分析，这种情况就类似于批次效应。我们生成一个随机的 $N \times 6$ 矩阵 $Y$ 来表示我们的研究对象，其中6表示不同的学科来表示 $N$ 表示 $N$ 个不同小孩。我们生成数据的方式就是，那些有一定相关性的学科。

1
2
3
4
5
6
7
8
9
10
library(MASS)
library(rafalib)
n <- 250
p <- 6
set.seed(1)
g <- mvrnorm(n,c(0,0),matrix(c(1,0.5,0.5,1),2,2))
Ystem <- g[,1] + matrix(rnorm(n*p/2,0,0.65),n,p/2)
Yhum <- g[,2] + matrix(rnorm(n*p/2,0,0.65),n,p/2)
Y <- cbind(Ystem,Yhum)
colnames(Y) <- c("Math","Science","CS","Eng","Hist","Classics")

样本相关性

我们要注意一下，这6个学科有着高度的相关性，如下所示：

1
round(cor(Y),2)

结果如下所示：

1
2
3
4
5
6
7
8
> round(cor(Y),2)
         Math Science   CS  Eng Hist Classics
Math     1.00    0.67 0.64 0.34 0.29     0.28
Science  0.67    1.00 0.65 0.29 0.29     0.26
CS       0.64    0.65 1.00 0.35 0.30     0.29
Eng      0.34    0.29 0.35 1.00 0.71     0.72
Hist     0.29    0.29 0.30 0.71 1.00     0.68
Classics 0.28    0.26 0.29 0.72 0.68     1.00

我们将使用图片来展示一下相关性，从下图我们可以发现，我们可以把学生分为STEM组与人文组（注：STEM与人文学科(humanities )其实就相当于中国的理科和文科，在这个案例中，STEM指的是科学(Science)，计算机(CS)和数学(Math)，人文学科指，古典学(Classics)，历史学(Hist)和英语语言文字(Eng)）。

在下图中，高度相关性的学科用红色表示，不相关的学科用白色表示，负相关的学科用蓝色表示，代码如下所示：

correlation_images,fig.cap
1
2
3
4
5
6
7
8
9
10
11
library(RColorBrewer)
mypar(1,2)
cols=colorRampPalette(rev(brewer.pal(11,"RdBu")))(100)
eps = matrix(rnorm(n*p),n,p)
par(mar = c(8.1, 8.1, 3.5, 2.1))
image(1:ncol(Y),1:ncol(Y),cor(Y)[,6:1],xaxt="n",yaxt="n",col=cols,xlab="",ylab="",zlim=c(-1,1),main="Actual Data")
axis(1,1:ncol(Y),colnames(Y),las=2)
axis(2,1:ncol(Y),rev(colnames(Y)),las=2)
image(1:ncol(Y),1:ncol(Y),cor(eps)[,6:1],xaxt="n",yaxt="n",col=cols,xlab="",ylab="",zlim=c(-1,1),main="Independet Data")
axis(1,1:ncol(Y),colnames(Y),las=2)
axis(2,1:ncol(Y),rev(colnames(Y)),las=2)

从上图中，我们可以看到以下信息：所有学科之间都有相关性，这表明学生成绩的背景有着某些隐藏因素（例如学术能力），正是这种能力导致了这些学科表现的相关性，因为在一个科目中获得高分的学生往往在其它科目也能获得高分。我们还能看到，这种相关性在STEM类学科之间更强，以及在人文学科之间更强。这就说明了，还有可能存在着另一种隐藏因素，这种因素决定了学生擅长STEM还是人文学科。现在我们使用一个统计模型来解释这种思路。

因子模型

基于前面的图形展示结果，我们假设这些现象背后有2个隐藏因素 $\mathbf{W}_1$ 和 $\mathbf{W}_2$ ，我们用这两个因素来解释所观察到的相关性结构，现在我们使用下面的公式来建模：

参数解释： $\alpha_{i,1}$ 表示学生 $i$ 的整体学术能力，$\alpha_{i,2}$ 表示学生 $i$ 的这种学术能力在STEM与人文学科之间的差异。现在我们的问题就是，如何估计 $W$ 和 $\alpha$ ？

因子分析与PCA

结果表明，在某些假设下，前两个主成分是对于 $W_{1}$ 和 $W_{2}$ 的最优估计，因此我们可以按以下方式对它们进行估计：

It turns out that under certain assumptions, the first two principal components are optimal estimates for $W_1$ and $W_2$. So we can estimate them like this:

1
2
3
4
s <- svd(Y)
What <- t(s$v[,1:2])
colnames(What)<-colnames(Y)
round(What,2)

结果如下所示：

1
2
3
4
5
6
7
> round(What,2)
      Math Science    CS  Eng Hist
[1,]  0.36    0.36  0.36 0.47 0.43
[2,] -0.44   -0.49 -0.42 0.34 0.34
     Classics
[1,]     0.45
[2,]     0.39

正如预期，第一个因子接近于一个常数，它有助于解释所有学科之间观察到的相关性，而第二个因子能解释STEM和人文学科之间的不表现。现在我们就以下模型中使用这些估计值：

我们现在就可以拟合模型，并注意一下，它能解释多大比例的变异：

1
2
fit = s$u[,1:2]%*% (s$d[1:2]*What)
var(as.vector(fit))/var(as.vector(Y))

结果如下所示：

1
2
3
> fit = s$u[,1:2]%*% (s$d[1:2]*What)
> var(as.vector(fit))/var(as.vector(Y))
[1] 0.7880933

我们在这里学习到的是，当我们有一些相关的单位或数据时，并不适合使用标准的线性模型。我们需要以某种试来解释观察到的结果。而因子分子则是实现这一目标的强大方法。

因子分析的常规用法

在高通量数据中，我们经常会看到相关结构。例如，我们在下图中展示了样本之间复杂的相关性。下面的这张图是按日期排序后，基因表达实验的数据展示结果：

gene_expression_correlations, fig.cap
1
2
3
4
5
6
7
8
9
10
library(Biobase)
library(GSE5859)
data(GSE5859)
n <- nrow(pData(e))
o <- order(pData(e)$date)
Y=exprs(e)[,o]
cors=cor(Y-rowMeans(Y))
cols=colorRampPalette(rev(brewer.pal(11,"RdBu")))(100)
mypar()
image(1:n,1:n,cors,col=cols,xlab="samples",ylab="samples",zlim=c(-1,1))

使用两个因子并不足以对观察到的相关结构进行建模。但是，我们可以使用一个更加普遍的因子模型：

我们可以使用PCA来估计$\mathbf{W}_1,\dots,\mathbf{W}_K$。但是，在选择 $k$ 方面并不容易，这就是我们要研究问题的核心。在下面的章节里，我们会介绍一下，如何使用探索性数据分析来辅助我们找到这个 $k$ 。

练习

P448

使用因子分析来对批次效应建模

继续使用前面的数据集：

1
2
library(GSE5859Subset)
data(GSE5859Subset)

我们之前展示过一些图形，那些图形中包括了性别效应与月份效应的基因子集，但是，现在我们使用一张图片来展示样本与样本之间的相关性（所有基因），并且展示出数据的复杂结构，如下所示：

correlation_image, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
library(rafalib)
library(RColorBrewer)
library(genefilter)
sex <- sampleInfo$group
batch <- factor(format(sampleInfo$date,"%m"))
chr <- geneAnnotation$CHR
tt<-rowttests(geneExpression,batch)
ind1 <- which(chr=="chrY") #real differences
ind2 <- setdiff(c(order(tt$dm)[1:25],order(-tt$dm)[1:25]),ind1)
set.seed(1)
ind0 <- setdiff(sample(seq(along=tt$dm),50),c(ind2,ind1))
geneindex<-c(ind2,ind0,ind1)
mat<-geneExpression[geneindex,]
mat <- mat -rowMeans(mat)
icolors <- colorRampPalette(rev(brewer.pal(11,"RdYlBu")))(100)
mypar(1,2)
image(t(mat),xaxt="n",yaxt="n",col=icolors)
y <- geneExpression - rowMeans(geneExpression)
image(1:ncol(y),1:ncol(y),cor(y),col=icolors,zlim=c(-1,1),
       xaxt="n",xlab="",yaxt="n",ylab="")
axis(2,1:ncol(y),sex,las=2)
axis(1,1:ncol(y),sex,las=2)

我们已经看到了即假设月份能够解释批次效应效应，并且使用线性模型对批次效应进行校正后的方法，这种方法表现很好。但是，这种方法仍然有改进的空间。这有可能是由于月份仅仅是背后一些隐藏因子的外面表现，或者是说是一些实际上导致结构或样本相关性的因子。

什么是一个批次

现在我画出每个样本的日期图，其中颜色表示月份，如下所示：

what_is_batch, fig.cap
1
2
3
4
times <-sampleInfo$date 
mypar(1,1)
o=order(times)
plot(times[o],pch=21,bg=as.numeric(batch)[o],ylab="Date")

我们注意到，每个月不止有一天。天这个时间单位也会产生影响吗？我们可以使用PCA和EDA来尝试回答一下这个问题。以下是按日期排序后的第一个主成分图：

PC1_versus_time, fig.cap
1
2
3
4
5
6
s <- svd(y)
mypar(1,1)
o<-order(times)
cols <- as.numeric( batch)
plot(s$v[o,1],pch=21,cex=1.25,bg=cols[o],ylab="First PC",xlab="Date order")
legend("topleft",c("Month 1","Month 2"),col=1:2,pch=16,box.lwd=0)

“天”似乎与第1主成分高度相关，它可以在很大程度上解释变异：

variance_explained, fig.cap
1
2
mypar(1,1)
plot(s$d^2/sum(s$d^2),ylab="% variance explained",xlab="Principal component")

进一步的研究显示，前6个左右的主成分(PC)似乎都是由日期驱动的：

PCs_stratified_by_time, fig.cap
1
2
3
4
5
mypar(3,4)
for(i in 1:12){
  days <- gsub("2005-","",times)  
  boxplot(split(s$v[,i],gsub("2005-","",days)))
}

如果我们从数据中将前6个主成分移除后，我们使用t检验来看一下结果是什么样子：

1
2
3
D <- s$d; D[1:4]<-0 #take out first 2
cleandat <- sweep(s$u,2,D,"*")%*%t(s$v)
res <-rowttests(cleandat,factor(sex))

事实上，这确实消除了批次效应，但是这似乎也消除了我们感兴趣的大部分生物学效应。事实上，没有基因的Q值再小于0.1了，如下所示：

pval_hist_and_volcano_after_removing_PCs, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
library(qvalue)
mypar(1,2)
hist(res$p.value[which(!chr%in%c("chrX","chrY") )],main="",ylim=c(0,1300))
plot(res$dm,-log10(res$p.value))
points(res$dm[which(chr=="chrX")],-log10(res$p.value[which(chr=="chrX")]),col=1,pch=16)
points(res$dm[which(chr=="chrY")],-log10(res$p.value[which(chr=="chrY")]),col=2,pch=16,
       xlab="Effect size",ylab="-log10(p-value)")
legend("bottomright",c("chrX","chrY"),col=1:2,pch=16)
qvals <- qvalue(res$p.value)$qvalue
index <- which(qvals<0.1)
cat("Total genes with q-value < 0.1: ",length(index),"\n",
    "Number of selected genes on chrY: ", sum(chr[index]=="chrY",na.rm=TRUE),"\n",
    "Number of selected genes on chrX: ", sum(chr[index]=="chrX",na.rm=TRUE),sep="")

这种情况似乎就是过度校正了，因为我们现在恢复了更少的Y染色体基因，并且p值的直方图中小p值的数目更小了，使得分布变得不再均匀。因此性别因素可能与前几个主成分(PC)有关，这就出现了”把婴儿和洗澡水一起倒掉“的情况。

SVA

前面我们刚提到了我们所遇到的问题，即我们遇到了过度校正，以及消除了与感兴趣结果相关的变异，解决这个问题的方法就是同时使用感兴趣变量的协变量以及那些被认为是批次的模型共同进行拟合。代理变量分析(SVA, Surrogate Variable Analysis)就是这样的一种方法。

SVA的基本思想就是先估计因子，但要注意不要包括感兴趣的结果。为此，SVA使用了一种交互式方法，其中的每一行都赋予一个权重，这个权重能够量化基因与替代变量（而不是感兴趣结果）唯一相关的概率。然后在SVD中使用这些权重，将更高的权重赋予给感兴趣的结果无关，且与批次有关的行。下面是这两次迭代的计算过程。三个图片显示的是，数据乘以权重（对于基因的子集）、权重，以及估计的第一因子，代码如下所示：

illustration_of_sva,fig.height
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
library(sva)
library(limma)
mod <- model.matrix(~sex)
cind <- order( as.Date(sampleInfo$date) )
dates <- gsub("2005-","",sampleInfo$date)
weights=rep(1,nrow(y))
par(mar = c(4.1, 2.1, 3.5, 2.1), 
    mgp = c(1.5, 0.5, 0))
layout(matrix(c(1:6),nrow=2,byrow=TRUE),widths=c(5,1.5,5))
for(b in 1:2){
  image(1:ncol(mat),1:nrow(mat),t(mat[,cind]*weights[geneindex]),xaxt="n",yaxt="n",col=icolors,xlab="",ylab="")
  axis(side=1,seq(along=dates),dates[cind],las=2)
  abline(v=12.5)
  
  svafit <- sva(y,mod,B=b,n.sv=5)
  weights = svafit$pprob.gam*(1-svafit$pprob.b)
  
  surrogate <- svd( y*weights)$v[,1]#Weighted SVD
  
  image(matrix(weights[geneindex],nrow=1),xaxt="n",yaxt="n",col=brewer.pal(9,"Blues"))
  plot(surrogate[cind],bg=sex[cind]+1,pch=21,xlab="",xaxt="n",ylab="Surrogate variable",ylim=c(-.5,.5),cex=1.5)
  axis(side=1,seq(along=dates),dates[cind],las=2)
  abline(v=12.5)
  text(1,0.5,"June")
  text(13.5,0.5,"Oct")
  legend("bottomright",c("0","1"),col=c(1,2),pch=16)
}

这个算法会多次迭代这个过程（由B参数控制），并返回替代变量的估计值，这就类似于因子分子中的隐藏因子。在R中，我们使用sva()函数来实现这个过程。在这个案例里，SVA会为我们选择代理值或因子的数量，如下所示：

1
2
3
4
5
library(limma)
svafit <- sva(geneExpression,mod)
svaX<-model.matrix(~sex+svafit$sv)
lmfit <- lmFit(geneExpression,svaX)
tt<- lmfit$coef[,2]*sqrt(lmfit$df.residual)/(2*lmfit$sigma)

结果如下所示：

1
2
3
> svafit <- sva(geneExpression,mod)
Number of significant surrogate variables is:  5 
Iteration (out of 5 ):1  2  3  4  5

与之前的方法相比，这种方法有了提升：

pval_hist_and_volcano_sva, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
14
res <- data.frame(dm= -lmfit$coef[,2],
                  p.value=2*(1-pt(abs(tt),lmfit$df.residual[1]) ) )
mypar(1,2)
hist(res$p.value[which(!chr%in%c("chrX","chrY") )],main="",ylim=c(0,1300))
plot(res$dm,-log10(res$p.value))
points(res$dm[which(chr=="chrX")],-log10(res$p.value[which(chr=="chrX")]),col=1,pch=16)
points(res$dm[which(chr=="chrY")],-log10(res$p.value[which(chr=="chrY")]),col=2,pch=16,
       xlab="Effect size",ylab="-log10(p-value)")
legend("bottomright",c("chrX","chrY"),col=1:2,pch=16)
qvals <- qvalue(res$p.value)$qvalue
index <- which(qvals<0.1)
cat("Total genes with q-value < 0.1: ",length(index),"\n",
    "Number of selected genes on chrY: ", sum(chr[index]=="chrY",na.rm=TRUE),"\n",
    "Number of selected genes on chrX: ", sum(chr[index]=="chrX",na.rm=TRUE),sep="")

结果如下所示：

1
2
3
4
5
6
> cat("Total genes with q-value < 0.1: ",length(index),"\n",
+     "Number of selected genes on chrY: ", sum(chr[index]=="chrY",na.rm=TRUE),"\n",
+     "Number of selected genes on chrX: ", sum(chr[index]=="chrX",na.rm=TRUE),sep="")
Total genes with q-value < 0.1: 14
Number of selected genes on chrY: 5
Number of selected genes on chrX: 8

现在我们通过图形来展示SVA实现了什么，下面是原始数据集分解为性别效应，代理变量和独立噪声的

为了可视化SVA实现了什么，下面是原始数据集通过算法分解为性别效果、代理变量和算法估计的独立噪声的可视化结果：

different_sources_of_var, fig.cap
1
2
3
4
5
6
7
8
9
10
11
Batch<- lmfit$coef[geneindex,3:7]%*%t(svaX[,3:7])
Signal<-lmfit$coef[geneindex,1:2]%*%t(svaX[,1:2])
error <- geneExpression[geneindex,]-Signal-Batch
##demean for plot
Signal <-Signal-rowMeans(Signal)
mat <- geneExpression[geneindex,]-rowMeans(geneExpression[geneindex,])
mypar(1,4,mar = c(2.75, 4.5, 2.6, 1.1))
image(t(mat),col=icolors,zlim=c(-5,5),xaxt="n",yaxt="n")
image(t(Signal),col=icolors,zlim=c(-5,5),xaxt="n",yaxt="n")
image(t(Batch),col=icolors,zlim=c(-5,5),xaxt="n",yaxt="n")
image(t(error),col=icolors,zlim=c(-5,5),xaxt="n",yaxt="n")

练习

P459